Modellering i python 🐍

Torodd F. Ottestad

🗒️ Plan

  • Kort intro til modellering
    • modellere lineær vekst
  • Modellering i ekte datasett
    • Eksponentiell vekst - BNP i USA (med litt teori)
    • Logistisk vekst - T.B.A…

Lineær vekst

\[f(x) = ax + b\]

Lage algoritme

Stigningstal: \[a = \frac{\Delta y}{\Delta x}\]

Multipliserer begge sider med \(\Delta x\)

\[\Delta y = \Delta x \cdot a\]

Lage algoritme

Prøv sjølv

Formuler ein algoritme (med ord) for korleis me kan bygga opp eit sett med \(x\) og \(y\)-verdiar ut frå samanhengen \[\Delta y = \Delta x \cdot a\]

  1. Bestemme startpunkt
  2. Finne neste punkt på grafen
    • ny \(x\) blir forrige \(x + \Delta x\)
    • ny \(y\) blir forrige \(y + \Delta y\)
  3. Gjenta til me har nådd stor nok \(x\)

Gjere om til python

Bruk modelleringsalgoritmen til å teikne grafen til \(f(x)=2x-1\) for \(x\in[0, 5]\) i python.

import matplotlib.pyplot as plt
import numpy as np

# tal frå oppgåveteksten
xstart = 0
xslutt = 5
ystart = -1
a = 2

# set steglengde lik 0.1
dx = 0.1

# lagar arrays til funksjonsverdiane
x_verdiar = np.array([xstart])
y_verdiar = np.array([ystart])

# lagar variablar x og y
x = xstart
y = ystart

# algoritmen
while x < xslutt:
    y = y + dx*a
    x = x + dx

    x_verdiar = np.append(x_verdiar, x)
    y_verdiar = np.append(y_verdiar, y)

# plottar resultatet
plt.plot(x_verdiar, y_verdiar)
plt.show()

Gjere om til python

BNP i USA

Last ned fila GDPUS_nsa.csv og legg ho i samme mappe som jupyter notatblokka de jobbar i.

Oppstart

Startar med å importera nødvendige bibliotek

import pandas as pd

import numpy as np
import matplotlib.pyplot as plt

Innlesing av data

df = pd.read_csv("GDPUS_nsa.csv")
df.DATE = pd.to_datetime(df.DATE)

df.head()

Innlesing av data

DATE NA000334Q
0 1947-01-01 58497.0
1 1947-04-01 60672.0
2 1947-07-01 62196.0
3 1947-10-01 68250.0
4 1948-01-01 64258.0

Utplukking

Nøyer oss med 1. januar kvart år:

df = df[::4]
df.head()
DATE NA000334Q
0 1947-01-01 58497.0
4 1948-01-01 64258.0
8 1949-01-01 66185.0
12 1950-01-01 68032.0
16 1951-01-01 81603.0

Kun årstal

Treng heller ikkje datoane no. Endrar òg indeks til årstal.

df.DATE = pd.DatetimeIndex(df.DATE).year
df.set_index("DATE", inplace=True)
df.head()
NA000334Q
DATE
1947 58497.0
1948 64258.0
1949 66185.0
1950 68032.0
1951 81603.0

Series

hentar kun ut bnp-verdiane. Med index.

bnp = df.NA000334Q
bnp.head()
DATE
1947    58497.0
1948    64258.0
1949    66185.0
1950    68032.0
1951    81603.0
Name: NA000334Q, dtype: float64

Plotting

bnp.plot()
plt.title("BNP i USA")
plt.show()

Modellering av vekst

Eksponentiell vekst

Eksponentialfunksjonen er gitt ved

\[f(x)=a\cdot e^{kx}\]

Dersom me deriverer den får me

\[f'(x) = k \cdot a \cdot e^{kx}\]

Og dermed \[f'(x) = k \cdot f(x)\]

Legg merke til

Eksponentiell vekst har me når vekstfarten er proporsjonal med funksjonen:

\[f'(x) = k \cdot f(x)\] \[y' = k \cdot y\]

For å modellera dette kan me bruka Eulers metode. Den tar utgangspunkt i at

\[f'(x)\approx \frac{\Delta y}{\Delta x}\]

Og dermed er \[ \Delta y \approx \Delta x \cdot f'(x) \]

Eksponentialfunksjonen

Prøver å bruka Eulers metode for å modellera \(f(x)=e^x\)

Koden

import matplotlib.pyplot as plt
import numpy as np

x_0 = 0                     # startverdi x
y_0 = np.exp(x_0)           # startverdi y

dx = 0.1                    # steglengde

x_verdiar = np.array([x_0])
y_verdiar = np.array([y_0])

for i in range(int(5/dx)):
    y_1 = y_0 + dx * y_0
    y_0 = y_1
    x_verdiar = np.append(x_verdiar, dx*(i+1))
    y_verdiar = np.append(y_verdiar, y_0)

Koden forts.

x_funk = np.linspace(0,5,100)
y_funk = np.exp(x_funk)

plt.plot(x_verdiar, y_verdiar, "r:", lw=5)
plt.plot(x_funk, y_funk, "k", lw=3)
plt.title(f"Modellert $e^x$ med $\Delta x$ = {dx}", fontsize=16)
plt.show()

Koden forts.

Tilbake til BNP

Eksponentiell vekst i datasettet

Eulers metode \[\Delta y \approx \Delta x \cdot f'(x)\]


Vekstfarten til \(f(x)=a\cdot e^{k\cdot x}\) \[ f'(x) = k\cdot f(x)\]

xstart = df.index[0]
xslutt = df.index[-1]
ystart = bnp.iloc[0]

dx = 0.01
k = 0.08

x_verd = []
y_verd = []

x = xstart
y = ystart

while x<=xslutt:
    x_verd.append(x)
    y_verd.append(y)
    y = y + dx*k*y
    x = x + dx
plt.plot(bnp, label="Datasettet")
plt.plot(x_verd, y_verd, label="Modellen")
plt.title(f"Modell for BNP i USA frå {xstart} til {xslutt}", fontsize = 18)
plt.legend()
plt.show()

Litt upresis \(k\). Vil prøva å estimera ein ny \(k\) ut frå datasettet.

Tar utgangspunkt i at \(f(x)=a\cdot e^{kx}\). Då er

\[f(0)=a \cdot e^{k\cdot 0} = a\]

og

\[f(1) = a \cdot e^{k\cdot 1}\]

Dermed blir

\[\frac{f(1)}{f(0)} =\frac{a\cdot e^k}{a} = e^k\]

Og vidare

\[\ln\left(\frac{f(1)}{f(0)}\right) = \ln(e^k) \Leftrightarrow k = \ln\left(\frac{f(1)}{f(0)}\right)\]

Dermed kan me legga inn utrekning av \(k\) rett før cella der me køyrer modellen:

ratio = bnp.div(bnp.shift(1))
snittratio = np.average(ratio.dropna())
forslag_k = np.log(snittratio)

print(f" k = {forslag_k:.4f}")
 k = 0.0621
xstart = df.index[0]
xslutt = df.index[-1]
ystart = bnp.iloc[0]

dx = 0.01
k = forslag_k

x_verd = []
y_verd = []

x = xstart
y = ystart

while x<=xslutt:
    x_verd.append(x)
    y_verd.append(y)
    y = y + dx*k*y
    x = x + dx

plt.figure(figsize = (16, 9))
plt.plot(bnp, label="Datasettet")
plt.plot(x_verd, y_verd, label="Modellen")
plt.title(f"Modell for BNP i USA frå {xstart} til {xslutt}", fontsize = 26)
plt.legend()
plt.show()

Regresjon

Samanliknar modellen me fekk med modellen me ville fått ved regresjon.

Regresjon med scipy

# eksponentiell modell (generell)
def f(x, a, b):
    return a*b**x

# gjer om x-verdiane til "år etter 1947"
x_data = bnp.index.values - 1947
y_data = bnp.values

# startgjett [x, y]
p0 = [1, 0]

# får både parameterane P og kovarians K som output
P, K = curve_fit(f, x_data, y_data, p0)

# hentar ut a og b frå P
a, b = P

# lagar y-verdiar ut frå funksjonen (med a og b frå regresjonen)
y_reg = f(x_data, a, b)

# plottar og pyntar
plt.figure(figsize=(10, 6))
plt.scatter(x_data, y_data, label="Datasettet", color="blue")
plt.plot(x_data, y_reg, label="Regresjon", color="red")
plt.legend()
plt.xlabel("x, år etter 1947")
plt.ylabel("BNP i USD")
plt.show()

Regresjon med scipy

Logistisk vekst

Definisjon

Veksten er logistisk når vekstfarten følger modellen gitt ved

\[f(x)=k\cdot f(x) \cdot \left(1 - \frac{f(x)}{B}\right)\]

❓ eksponentiell ➡️ logistisk modell
Kva endringar må me gjera i koden for modellen til eksponentiell vekst for å kunna bruka han for å sjå på logistisk vekst?

Oljeproduksjon i Noreg

Ser på eit nytt datasett: oljeproduksjon-norge.csv

Kikkar på fila, den ser sånn ut i starten:

import matplotlib.pyplot as plt
import pandas as pd 
import numpy as np

# les inn datasettet
oljedata = pd.read_csv("oljeproduksjon-norge.csv", sep=";", decimal=".", comment="#")

# lagar eigne arrays for ti dog produksjon
tid = np.array(oljedata["År"])
produksjon = np.array(oljedata["Produksjon"])

# plottar datapunkta
plt.plot(tid, produksjon, "r.")
plt.show()

Logistisk vekst i datasettet

Ser av figuren at me kan anta bæreevne på ca. 250, dvs \(B=250\).

import matplotlib.pyplot as plt
import pandas as pd 
import numpy as np

# les inn datasettet
oljedata = pd.read_csv("oljeproduksjon-norge.csv", sep=";", decimal=".", comment="#")

# lagar eigne arrays for ti dog produksjon
tid = np.array(oljedata["År"])
produksjon = np.array(oljedata["Produksjon"])

# plottar datapunkta
plt.plot(tid, produksjon, "r.")

# Bygger modellen
xstart = tid[0]
xslutt = tid[-1]
ystart = produksjon[0]

dx = 0.01
k = 0.5 # TIPPING. prøv deg fram til grafen får rett form
B = 250 # LEST ut frå datapunkta

x_verdiar = []
y_verdiar = []

x = xstart
y = ystart

while x <= xslutt:
    x_verdiar.append(x)
    y_verdiar.append(y)
    y = y + dx * k * y * (1 - y/B) 
    x = x + dx

# teiknar modellen sin graf
plt.plot(x_verdiar, y_verdiar)
plt.show()

Kan tilpassa formen med å endra k=...
Kan tilpassa posisjon i x-retning med å endra xstart = tid[0]